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Direct Multi-grid Methods for Linear Systems 
with Harmonic Ahasing Patterns 

Pablo Navarrete Michelini 



Abstract 

Multi-level numerical methods that obtain the exact solution of a linear system are presented. The 
methods are devised by combining ideas from the full multi-grid algorithm and perfect reconstruction 
filters. The problem is stated as whether a direct solver is possible in a full multi-grid scheme by 
avoiding smoothing iterations and using different coarse grids at each step. The coarse grids must form 
a partition of the fine grid and thus establishes a strong connection with domain decomposition methods. 
An important analogy is established between the conditions for direct solution in multi-grid solvers and 
perfect reconstruction in filter banks. Furthermore, simple solutions of these conditions for direct multi- 
grid solvers are found by using mirror filters. As a result, different configurations of direct multi-grid 
solvers are obtained and studied. 

Index Terms 

multigrid, perfect reconstruction filter, domain decomposition, direct solver, aliasing. 

I. Introduction 

This study focuses on the problem of solving the linear system of equations 

Au = f (1) 

over the field of complex numbers. The problem is restricted to the case when the number of equations, 
n, is the same as the number of unknowns, and the number n is even (in some cases a power of 2). 
The system matrix A G C"^" is sparse and it will be assumed to be invertible, with special attention to 
ill-conditioned cases. The problem becomes challenging when n scales to large numbers (e.g. thousands 
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of unknowns). This situation arises frequently in scientific and engineering computations, most notably 
in the solution of PDEs [1], [2] and other areas Uke simulation of stochastic models [3] and the solution 
of optimization problems [4]. 

A vast amount of numerical methods exist to solve this problem efficiently. They vary from direct 
(or exact) solvers that compute the exact solution, u, to iterative solvers that compute a sequence of 
approximations that converges to the exact solution, Vk —>■ u. Here, convergence must be defined to 
minimize some norm of the approximation error, = u — v^- Direct solvers are often based on some 
type of matrix factorization, being the most popular those based on LU decomposition [5], [6]. Among 
the iterative solvers, the most common are: stationary iterative methods (e.g. Gauss-Seidel, Jacobi and 
Richardson iterations), Krylov subspace methods (e.g. conjugate gradients, GMRES and BiCG) and 
multi-level methods (e.g. multi-grid and domain decomposition) [1], [2]. 

In this study, multi-level numerical methods working as direct solvers are obtained. In the same category 
there are other direct multi-level solvers like: total reduction methods [7], [8], partial (cyclic) reduction 
methods [9] and LU factorization of non-standard forms [10]. Besides their structural differences, 
each method works under certain limitations. Total reduction and partial (cyclic) reduction methods 
are specifically designed for Poisson's equation, and LU factorization of non-standard forms works for 
elliptic problems. In this study, the limitations are not described in terms of categories of PDEs but in 
terms of two additional properties on the system. These are: A has to be diagonalizable, and ignoring 
some of the unknowns should produce a specific aliasing pattern. 

First, by assuming that the system matrix is diagonalizable, we have the eigendecomposition 

A = WAV^ , (2) 

where the columns of W form the set of right eigenvectors, A is a diagonal matrix with the eigenvalues 
of A, and the columns of V form the set of left eigenvectors. Here, the right and left eigenvectors form 
a biorthogonal basis so that V^W = / and they do not need to be equal, which allows the case of a 
non-symmetric system matrix. This restriction limits the applications to non-defective problems. 

The second restriction is on the effects of down-sampling the eigenvectors of the system. This operation 
drops a number of components when applied to vectors and keeps the remaining components untouched. 
An example is shown in Fig. [T] where one every two samples of harmonic functions are dropped. By 
down-sampling the eigenvectors of the system their linear independence is lost, because the dimension 
of the down-sampled space is less than the number of eigenvectors. This is a general description of the 
phenomenon of aliasing. Here, it will be assumed that there is a subset of n/2 components defining 
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(e) ws and down-sampled ws (f) W7 and down-sampled W7 (g) wg and down-sampled we (h) ws and down-sampled ws 

Fig. 1. Example of harmonic aliasing pattern in a harmonic sine basis. The set of vectors Wj £ R*, j = 1, ... ,9, with 
components (lUj)i = |sin(^), i = 1,...,9 forms an orthonormal basis of R*. A down-sampling operation drops the 
components with even i. The resulting down-sampled eigenvectors are not linearly independent. The down-sampling of Wk is 
equal to the down-sampling of Wg-k with fc = 1, . . . , 4. 



a down-sampling operation that makes each down-sampled eigenvector equal (up to sign) to only one 
of the other down-sampled eigenvectors (see Fig. [T|). This specific pattern is called harmonic aliasing 
pattern ril] and is found in harmonic functions, which are eigenvectors of linear space invariant (LSI) 
system^j, and other basis including at least: Hadamard matrices and eigenvectors of coupled systems of 
equations [11], [12]. 

The generality of a system with harmonic aliasing patterns is not known at its largest extent. They were 
introduced in [11] in order to extend the strong convergence analysis of multi-grid algorithms based on 
local Fourier analysis (LFA). This analysis was introduced by Achi Brandt in the late 70's and remains 
as the main rigorous tool for the design of multi-grid methods [13], [14]. LFA is based on Fourier 
analysis and is thus restricted to LSI systems, which makes it difficult to use it in many applications. 
The extended convergence analysis in [11] has not overcome this problem drastically. Nevertheless, in 
this study its algebraic framework will allow not only the analysis of traditional multi-grid methods but 
also the introduction of new direct solvers with connections to other important numerical methods. 

The contributions of this paper are thus both practical and conceptual. From the practical side, multi- 
grid algorithms have been successfully applied in many practical problems but the theory behind does not 
reach the same level. The most common implementation is algebraic multi-grid (AMG) which obtains a 
multi-grid configuration based on heuristics [15]. The numerical methods obtained in this paper represent 
a step forward on what the theory can achieve. This is, under the assumptions just mentioned, a completely 
algebraic configuration is found that solves the problem exactly with no use of heuristics. The algorithms 



'in numerical analysis a different terminology is used. The stencil of an unknown is defined as a geometric arrangements of 
the non-zero coefficients in the correspondent row in A, centered at the diagonal element. This is equivalent to the concept of 
impulse response in signal processing and an LSI system is equivalent to a system with constant stencil coefficients. 
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are computationally efficient and adaptable to the computational resources (single-core or multi-core). 

The conceptual contribution of this paper is the introduction of numerical methods in a setting that 
establishes a clear connection between multi-rate systems and different multi-level numerical methods. 
The analysis exploits the structural similarities of the full multi-grid algorithm [16] and problems of signal 
reconstruction in multi-rate systems [17]. On one hand, both classical and extended convergence analysis 
for multi-grid solvers have shown the importance of aliasing phenomena to explain how the algorithm 
converges [11], [15]. Here, aliasing appears as an additional source of error that has to be controlled 
by the algorithm. On the other hand, the problem of signal reconstruction in multi-rate systems follows 
a different approach. A signal is decomposed in different coarse levels and the question is posed as 
whether the original signal can be reconstructed from these different pieces of information. Here, each 
coarse level component has aliasing effects and perfect reconstruction is possible because these effects 
cancel each other. Given the structural similarities between multi-grid methods and multi-rate system, 
the central question is whether aliasing cancellation can be used in multi-grid to obtain direct solvers, 
where the analogue of "perfect reconstruction" is "exact solution." 

The main goal, under the assumptions just mentioned, is to modify a full multi-grid algorithm and 
obtain a direct multi-grid solver in direct analogy with the problem of perfect reconstruction filters. In 
order to establish this analogy the problem of perfect reconstruction needs to be generalized for systems 
with harmonic aliasing patterns, which allows to configure perfect reconstruction filters that are not 
necessarily LSI systems. On the other hand, the full multi-grid algorithm also requires modifications. 
Multiple coarse grids are needed in order to keep all the information from the original problem in coarse 
levels. A partition of the complete set of unknowns defines several coarse levels and represents a particular 
type of domain decomposition [18]. 

Similar approaches can be found in the literature. The use of multiple coarse grids in multi-grid has 
been introduced by Frederickson and McBryan in [19] and Hackbusch in [20], [21]. It is known that 
aliasing cancellation helps these methods to converge fast [22]. Nevertheless, none of them work as direct 
solvers and they still use smoothing iterations. Total reduction methods are the closest in structure to the 
direct solvers obtained and can be seen as a more restrictive version of one of the algorithms presented. 

In section JIJ full two-grid algorithms and multi-rate systems are reviewed. In section JIIJ harmonic 
aliasing patterns are introduced. In section |IVl perfect reconstruction filters are studied in the context of 
harmonic aliasing patterns. In section |Vl the convergence analysis of two-grid methods based on grid 
partitions is studied. In section |Vll the problem of finding direct two-grid solvers is stated and solved. 
In section IVIIl the multi-grid case is considered. Finally, in section IVIIII some examples are presented. 
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(a) Full two-grid algorithm (b) Two-channel multi-rate system 



Fig. 2. The numerical methods introduced in this paper are based on two systems. First, the full two-grid algorithm in Fig. [2a] 
This is an iterative algorithm to obtain an approximate solution of Au = /. The dotted line separates vectors from the fine and 
coarse grid domains. The interpolation (restriction) operation is applied to vectors crossing the dotted line from below (above). 
Second, the two-channel multi-rate system in Fig.|2b] Here, each box represents an LSI system and inside the stationary impulse 
response is shown. The circles with "| 2" represent down-sampling operators that drop one every two samples. Similarly, circles 
with "I 2" represent up-sampling operators that insert one zero every two samples. 



II. Preliminaries 

A. Full two— grid algorithm 

The full two-grid algorithm is an iterative solver used to obtain approximate solutions of ([T]). In order 
to simplify the problem, the algorithm uses the concept of grids and coarse grids. Whatever the nature 
of the problem is, the system ([T]) can always be associated with a graph in which the unknowns of the 
system are the nodes of the graph. The nodes of the graph are associated with a set of labels O which 
is called the fine grid. A coarse grid, Q, is a proper subset of the fine grid; i.e., 1^ C 1^. 

The so-called inter-grid operators are defined as any linear transformation between scalar fields on 
and 0. That is 7/ G £,nx\n\ ^ (j--|r2|xn^ where Ij is the interpolation operator and Jr is 

the restriction operator. In addition to these operations, and following the standard of most multi-grid 
applications, a coarse system matrix is defined following the Galerkin condition [16] 

A'^IrAIi. (3) 

Kfiill two-grid algorithm solves (d) by using the system shown in Fig. [2a| Here, there are three steps 
involved. First, the two boxes perform fixed numbers of smoothing or stationary iterative iterations. 
These iterations -typically Gauss-Seidel, Jacobi, Richardson, etc.- are known to obtain good local 
approximations of the solution [16]. This means that high-frequency components of the approximation 
error, = u — Vk, are efficiently reduced. In each smoothing iteration the approximation error evolves 
as Cfc+i = Sck. The matrix S is thus called the smoothing operator. For better understanding of this 
step is convenient to assume that 5 is a filter, or Fourier multiplier. If this is the case then S has same 
eigenvectors as A and it can be decomposed as 5" = WEV^ . The matrix of eigenvalues, S, contains the 
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s = w 



V 



H 



(4) 



frequency response, or symbols, of the smoothing operator. The smoothing effect on the approximation 
error is seen in the frequency domain as a damping effect concentrated on high-frequency eigenvectors. 
In other words, the eigen-decomposition of 5 can be written in block-form as 

where and Y^h correspond to low- and high-frequency eigenvalues close to 1 and 0, respectively. 

The remaining steps in Fig. |2a] make use of the coarse grid. First, the so-called nested iteration step 
takes / and computes an initial approximation vq. This step solves a coarse grid equation using a restricted 
version of the source vector, /. And second, the so-called correction scheme improves the approximation 
vi to obtain V2- This step computes an approximation of tk from the error equation Ack = Vk, where 
Vk = f — Avk is the residual vector taking the role of the source vector. A coarse grid equation is solved 
using a restricted version of the source vector, fk- Once the approximation is obtained, it is added to 
correct the current approximation. 

It is observed that the approximation of obtained in the coarse grid equations effectively represents its 
low-frequency components and fails to represents its high-frequency components. This is a consequence 
of the computations in coarse grids where the big picture (low-frequencies) of the solution is clear but 
details (high-frequencies) are lost. 

The approximation error in the coarse grid steps evolves as 

eo = Ku and e^+i = Kek , (5) 
for nested iterations and the correction scheme, respectively. Here, the matrix 



K = I- IiA-^IrA 



(6) 



determines the evolution of the approximation error and is called the coarse grid correction matrix [23]. 

The reduction of low-frequency components of the error in two-grid steps suggests that K is a high- 
pass filter. If this would be the case then there would be a decomposition K = WTV^ where F is a 
diagonal matrix. One would expect a frequency response of the filter with the shape shown in Fig. [3al 
Unfortunately, this is not the case. K is not a filter because of aliasing effects. 

Although K is not technically a filter (or Fourier multiplier), under the assumption of harmonic aliasing 
patterns a decomposition K = WFV^ exists where f is sparse (but not diagonal). In the forthcoming 
sections it will be shown that 



K = W 



V 



H 



(7) 
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(a) Two-grid filtering effect (b) Two-grid aliasing effect 

Fig. 3. Representation of the frequency effects of a coarse grid correction operator in a two-grid algorithm. The operator is 
not a filter and thus is not completely represented by filtering effects. The operator is completely represented by filtering and 
aliasing effects. Tl~>l indicates the region of low-frequencies where the symbols are close to 0. f indicates the region 
of high-frequencies where the symbols are close to 1. Tl^h and Th^l multiply low- and high-frequencies and take them 
into high- and low-frequencies, respectively. 



where Tl^l, ^h^l, ^l^h and Th~>h are all diagonal matrices if one assumes harmonic aliasing 
patterns. The filtering effect is contained in T^^l and Th^h, and are expected to be as shown in Fig. 
l3al As opposed to a proper filter, this figure does not tell everything about the coarse grid correction 
matrix. A second graphic, shown in Fig. [3bl must show the aliasing effect from Th^l and Tl^h. 

B. Two-channel multi-rate systems 

In multi-rate systems one is interested to decompose a discrete signal in different components at 
lower sampling rates [17], [24]. A two-channel multi-rate system is shown in Fig. [20 Two restriction 
operations are performed by filtering and down-sampling. These operations split the original signal, s, 
into two signals, sq and si, each with half of the original samples. The original signal is recovered by 
summing two interpolation operations performed by inserting zeros (up-sampling) and then filtering. 

Here, the boxes represent LSI systems which have stationary impulse responses and are filters with 
respect to a harmonic basis. Their frequency responses are given by the Fourier transforms of their 
impulse responses: Hq{uj), Hi{uj), Go{lo) and Gi{uj). In practical applications the problem is whether 
the support of the frequency responses can overlap and still be able to recover the original signal. This 
is possible when the filters fulfiU the following conditions by Vetterli [25] 

Go{io)Ho{u;) + Gi{Lj)Hi{io) = 2, (8) 
Goiio)Ho{u; - it) + Giiio)Hi{u; - tt) = . (9) 

Here, condition Q causes the aliasing effects from different channels to cancel each other and dD causes 
the final sum to be equal to the original signal. The more general result with an arbitrary number of 
decompositions is due to Vaidyanathan [26]. 
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(a) Fine grid witii nodes colored red and black. (b) Red and black coarse grids 

Fig. 4. Red-black partition of a grid in ID. The fine grid f2 is partitioned into a red coarse grid Q and a black coarse grid Q. 



III. Red-Black Harmonic Aliasing 

A red coarse grid, Cl, and a black coarse grid, Cl, are defined by a partition of the fine grid. This is, 

= 17 + il, where the sum denotes the disjoint union of two sets. An example is shown in Fig. HI The 
motivation of the red-black partition is to keep track of all the fine grid nodes in coarser grids. In this 
way the partition represents a particular type of domain decomposition [18]. 

The selection of nodes to the red and black partition will be represented by down-sampling operators 
according to the following definition. 

Definition 1 (Down/Up— sampling matrices): The red and black down— sampling matrices are defined 
as D e {0, l}|f^|x" and D G {0, Ijl^^lx" such that 



, f I 1 if node i G 17 is the i*^ red node in 17 
{D\,j { (10) 
otherwise 



and 



, 1 if node 7 G 17 is the i*'^ black node in 17 

my={ , (11) 

otherwise 

respectively. Similarly, the red and black up-sampling matrices are defined as [/ G {0,1}"^ 1^1 and 
U G {0, l}"x|f^l such that U = and U = , respectively. 

The down-sampling matrix of a certain color represents a linear transformation that takes a fine grid 
vector and drops all the values that correspond to a node of different color. The up-sampling matrix 
takes a coarse grid vector and inserts zeros at the new nodes (of different color) in the fine grid. From 
this interpretation, a set of basic properties follows 

DtJ = I and DU = I , (12) 
DU = and DJJ = , and (13) 
tJD + UD = I . (14) 
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(a) Low-frequency harmonic function wl- (b) High-frequency harmonic function wh- 




(c) DwL (d) DwL (e) Dwh (f) Dwh 

Fig. 5. Difference between red and black harmonic aliasing patterns. The red down-sampling of tol and wh are equal, 
Dwl = Dwh- Whereas the black down-sampling of wl is the negative of the black down-sampling of wh, Dwl ~ —Dwh- 



In [11] a so-called harmonic aliasing pattern was defined only for the red grid. The following definition 
considers both red and black grids. This addition has important consequences in the results to come. 

Definition 2 (Red and Black Harmonic Aliasing Patterns): A matrix M € C"^" is said to have red 
and black harmonic aliasing patterns if it is diagonalizable, with biorthogonal eigenvectors W and V, 
and there exists a red-black partition which divides the domain into two halves, with down-sampling 

_ n ^ n 

matrices D G {0, 1}2 and D G {0, 1}2 such that 



V^UDW = N and V"UDW = N , 



(15) 



respectively. Here, A'^ and N are the red and black harmonic aliasing patterns defined, respectively, as 

(16) 



I I 
I I 



and N = - 



I -I 
-I I 



In this definition the red and black harmonic aliasing patterns appear as independent properties. The 
following statement shows how these definitions are equivalent. 

Proposition 1: A matrix with red harmonic aliasing pattern has black harmonic aliasing pattern, and 
vice versa. Therefore, a matrix with these properties is said to have a red-black harmonic aliasing pattern. 
Proof: Taking ([141 ). pre-multiplied by and post-multiplied by W gives 



V^UDW + V^UDW = I . 



(17) 



Thus, if V^UDW = N then V^UDW = N, and vice versa. ■ 
The definition of red-black harmonic aliasing pattern is convenient for algebraic manipulation but the 
connection with the common concept of aliasing is not clear yet. In the following theorem an alternative 
and equivalent definition is given which makes this connection explicit. 
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Theorem 1 (Navarrete and Coyle): A matrix M £ C"^" with red-black harmonic aUasing pattern is 
equivalent to have a diagonalizable matrix, with biorthogonal eigenvectors W and V, for which there exists 
a red-black partition dividing the domain into two halves, with down-sampling matrices D G {0, 1} 2 " 

~ n 

and D G {0, 1} 2 and such that there is an ordering of the eigenvectors for which the partitions 
W = [WlWh] and V = [VlVh] fulfill the conditions 

DWl = DWh , DVl = DVh , (18) 

DWl = -DWh and DVl = -DVh- (19) 

The proof of this theorem is partially contained in [11] were only the red coarse grid was considered. 
The result for the black coarse grid is shown in Appendix [A] 

The sign difference between the red and black harmonic aliasing patterns is explained in Fig. [5] and 
it represents the fact that harmonic basis vectors are composed of two envelopes which, intermixed with 
the same sign form a low frequency and, intermixed with opposed signs form a high-frequency. 

The definition of red-black harmonic aliasing patterns for a given matrix does not involve its eigen- 
values. But, in the algebra derived from the partition of eigenvectors in Theorem [T] it will be necessary 
to specify which eigenvalues are associated with each partition. The following remark introduces the 
notation to make this distinction clear. 

Remark 1: For a matrix M G C"^" with red-black harmonic aliasing patterns and eigen-decomposition 
M = WEV^, the partition of eigenvectors, W = [WlWh] and V = [VlVh], induces a partition 
of eigenvalues such that E = Eh ' where El and Eh are the diagonal matrices of eigenvalues 
associated with L and H eigenvectors, respectively. 

IV. Perfect Reconstruction Filters for Systems with Harmonic Aliasing Patterns 

In the context of finite discrete systems we want to extend the idea of perfect reconstruction filters for 
systems with harmonic aliasing patterns, which are more general than LSI systems. First, we define an 
object that will extend the operation of quadrature and conjugate mirror filters [27]-[29]. 

Definition 3 (Mirror Matrix): The mirror of a matrix M G C"^" with respect to a red-black partition 
represented by down/up-sampling matrices D, U = D^, and D, U = W , is defined as 

M* = itJD - ub)M{UD - UD) . (20) 

Matrices J7Z) G {0, and UD G {0, 1}"^" are diagonal with one's whenever z = j is a red or 

black node, respectively, and zero otherwise. Therefore, tJD — UD is a diagonal matrix that takes the 
value 1 when i = j is a red node, and takes the value —1 when i = j is a. black node. 
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Fig. 6. Two-channel finite multi-rate system. The system takes a signal s G C", decomposes it into s G C"^^ and s £ C"^'^, 
and recovers t £ C". A perfect reconstruction system is such that t — s. Here, Fa and Fr are restriction filters; Fi and Fj 
are interpolation filters; D and D represent down-sampling operations; and, U and U represent up-sampling operations. 



For an LSI system with stationary impulse response h[n\, the mirror operator with respect to a uniform 
down-sampling by factor of 2 (where red nodes are even nodes) produces an impulse response h* [n] = 
(— l)"/i[n]. This filter has symbols H*{lu) = H{lo — vr) and thus swaps low and high frequencies. The 
following proposition generalizes this result to systems with red-black harmonic aliasing patterns. 

Proposition 2: The mirror of a matrix M e ^"X" ^vith red-black harmonic aliasing patterns and 



eigen-decomposition M = W 



, is a filter with eigen-decomposition 



M" = W 



Eh 



El 



V 



H 



(21) 



Proof: Using the eigen-decomposition of M and the definition of red and black harmonic aliasing 
patterns, the result is obtained as follows 

M* = WV^{UD - UD)W E V"{UD - UD)WV^ 



El 
Eh 



(22) 



W 



V 



H 



These results are all we need to proceed into the results of the next sections. Now, in order to show the 
analogy between perfect reconstruction filters and direct multi-grid solvers we have to look back to the 
problem of perfect reconstruction and see how the main results look for systems with harmonic aliasing 
patterns. First, we need to restrict the problem to discrete signals in C". The two-channel multi-rate 
system in Fig. [6] is then the finite and discrete version of the system in Fig. [20 The problem is to 
find interpolation and restriction matrices Fj, Fr, Fj and in C"^" that are filters with respect to a 
biorthogonal basis W and V, and such that we have perfect reconstruction; i.e., t = s. 

The following theorem restates Vetterli's conditions ^ and ^ for this new context. 
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Theorem 2 (Vetterli): Let Fj, Fr, Fj and Fr be filters with respect to a biorthogonal basis W and 
V in C"^". Let their matrices of eigenvalues be H/, IIr, 11/ and II/j, respectively, and the eigenvectors 
have red-black harmonic aliasing patterns with respect to the down-sampling matrices D and D. Then, 
the multi-rate system in Fig. [6] has the perfect reconstruction property, t = s, if and only if 

tlitlR + UiUr = I , (23) 

fii,Lf^R,H - f^I,Lf^R,H = , and (24) 

til,HfiR,L - fll,HflR,L = , (25) 

where the L and H subindexes follow the notation introduced in Remark [T] 
Proof: From Fig. [6l perfect reconstruction is obtained if and only if 

FiUDFr + FiUDFr = I . (26) 

Pre-multiplying by , post-multiplying by W and using the definitions in (fTSl) gives 

^I,H^R,L—^I,H^R,L ^I,H^R,H-\-Tll^H^R,H 

The blocks in the diagonal are equivalent to (l23l) and the off-diagonals give (l24l) and (l25l) . ■ 
Next, we are interested in a particular solution of these conditions using mirror filters. The following 

corollary generalizes the solutions using quadrature and conjugate mirror filters (QMF and CMF) for 

perfect reconstruction in a two-channel multi-rate system with respect to a biorthogonal basis with 

red-black harmonic aliasing patterns. 

Corollary 1 ( Croisier et al. / Smith and Barnwell III / Mintzer): Perfect reconstruction conditions (1231) . 

(I24l l and ^ are fulfilled if Fj is such that 

n|,L + nf/, = I, and (28) 

Fr = Fi , Fi = F/ , and Fr = F/ . (29) 

Proof: Using Theorem |2j if ( [29l ) is assumed then the conditions (l24b and (|25] ) are fulfilled. The 
mirror condition (|28] ) is necessary to fulfill (|23] ). ■ 
Quadrature mirror filters were introduced as a solution of the perfect reconstruction problem in [27]. In 
its original formulation, only a Haar filter gives a sparse quadrature mirror filter [30]. Later, this problem 
was solved by the introduction of conjugate mirror filters [28], [29], and later generaUzed for biorthogonal 
and multi-channel paraunitary systems [25], [26]. The problem to obtain perfect reconstruction FIR filters 
has to do with the delays between the filters. From a signal processing perspective the red-black partition 



"7 0" 
.0 I- 



(27) 
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corresponds to a polyphase decomposition [17], [24], [31] that introduces the delays in the down/up- 
sampling operations and allows to condense the different solutions for perfect reconstruction. 

V. Two-grid Convergence 

Since the configuration of a two-grid algorithm involves many parameters and assumptions, it is 
convenient to introduce a single terminology that refers to this configuration. 

Definition 4 (Red— Black Harmonic Two— grid Configuration): A red— black harmonic two— grid config- 
uration is a set of matrices depending on a system matrix A G C"^" with red-black harmonic aliasing 
patterns and biorthogonal eigenvectors W and V. The red-black partition of nodes is represented by 
down-sampling operators D and D. The configuration is completed with the inter-grid filters Fr, Fj, 
Fji and Fj, all in C"^". Then, a series of matrices associated with the configuration are defined. First, 
the interpolation and restriction operators, and the two coarse system matrices given by the Galerkin 
condition, are defined as 



Ir = DFr , 
A = IrAIi 



h = FiU 
and 



Ir = DFr , 
A = IrAIi . 



Ii = FiU , (30) 
(31) 



The system matrix and the inter-grid filters have eigen-decompositions 



A = W AV" , Fr = WUrV 



H 



Fi = WUiV 



H 



Fr = WUrV" and Fi = WUiV" 



and the partition of eigenvectors, W = [WlWh] and V = [VlVh], leads to the partitions of eigenvalues 



A 



Hi? 



Hr,. 



and 11/ 



fll,H 



Finally, the two coarse grid correction matrices are defined as 

K = I- IjA-^IrA and 
K = I- IiA-'IrA . 



(32) 
(33) 



Now, based on these definitions, the goal is to obtain an eigen-decomposition of the coarse grid 
correction matrices. The following lemma takes the first step by giving useful expressions for the inverse 
of the coarse system matrices in (|32] ) and (1331) . 
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Lemma 1 (Navarrete and Coyle): In a red-black harmonic two-grid configuration the inverses of the 
red and black coarse grid matrices are given by 



A-^ = 4 {DWl)A-\DVl)" and 



A-^ = 4 iDWL)A-\DVLf 



(34) 
(35) 



respectively, with 



A = Un,LALni,L + TiR,HAHni,H and (36) 

A = Ur^lAlUj^l + flR,H^Hf^I,H ■ (37) 

The proof of this lemma is partially contained in [11] were only the red coarse grid was considered. 
The proof for the black grid is shown in Appendix iBl 

The result in Lemma [T] reflects the structure of a system with harmonic aliasing patterns. The eigen- 
vectors of coarse system matrices are given by a linear-independent subset of the down-sampling 
eigenvectors of the system matrix. The coarse eigenvalues are expressed as sums of low- and high- 
frequency eigenvalues as a result of aliasing effects. 

Finally, the following theorem gives the eigen-decompositions of coarse grid correction matrices. 

Theorem 3 (Navarrete and Coyle): In a red-black harmonic two-grid configuration the red and black 
coarse grid correction matrices (l32l) and (l33l) can be decomposed as 



K = W 



and 



(38) 



K = W 



V 



H 



(39) 



with 



and 



Tl^l = i- n/,LA-^nR,iAL 



Tl^h = IIi^hA-^IIr,l^l 



Th^H = I- Ill,H^-^t^R,H^^H 

Th^L = fll,L^-^fiR,Hf^H , 
Th^H = I- fll,H^-^flR,HAH . 



The proof of this theorem is partially contained in [11] were only the red coarse grid was considered. 
The result for the black coarse grid is shown in Appendix O 

The results for the red and black coarse grids carry the difference in sign from the definition of 
harmonic aliasing patterns, which can be seen in the cross-modal symbols (H L and L H). 
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VI. Direct Two-grid Methods 

In this section the full two-grid algorithm shown in Fig. |2a] will be modified to obtain direct two-grid 
solvers for systems with harmonic aliasing patterns. The motivation is to eliminate smoothing iterations 
in the full two-grid scheme and base the algorithm purely on nested iterations and/or correction schemes. 

A mere elimination of smoothing iterations in Fig. |2a] would make it impossible for the algorithm to 
converge since partial information in a single coarse grid is not enough to get all the information from 
the fine grid. The algebra is very clear on this point because, based on the Galerkin condition, a coarse 
grid correction matrix is idempotent (or projection matrix). This is, 

= K and K'^ = K , (40) 

which means that several iterations of two-grid steps with a single coarse grid do nothing more than 
a single iteration. On the other hand, the red-black partition keeps all the information from the fine 
grid. Therefore, combining red and black coarse grids at different steps of the algorithm has a chance to 
converge depending on the configuration of the algorithm. 

A. Multiplicative Approach 

Two modifications of the full two-grid algorithm from Fig.|2a]are considered in Fig.jTa] First, smoothing 
iterations are removed. And second, red and black coarse grids are considered in the nested iteration and 
correction scheme steps, respectively. 

If this algorithm works as a direct solver then it means that the problem is being factorized into coarse 
grid sub-problems. The following proposition shows this in terms of a factorization of A~^. 

Proposition 3: The system shown in Fig. ITal works as a direct solver; i.e., v = u, if and only if the 
following decomposition applies on the inverse of the system matrix: 

A-^ = IiA-^Ir + IiA-^Ir - IjA~^ {IrAIi)A-^Ir . (41) 

Proof: Using ([5]l in Fig. |7a] gives 

e = KKu . (42) 

Then, an exact solution is obtained if and only if KK = 0. Using (|32l ) and (|33] ) gives (|4TI ). ■ 
In this matrix factorization, the first two terms indicate the components of the inverse that come from 
the red and coarse grids independently. The third term indicates the dependence between the two solutions. 
In fact, the correction scheme in Fig. |7a] works on top of the solution given by nested iteration and thus 
mixes the solutions from both coarse grids. 
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This factorization is known in domain decompositions as the multiplicative Schwartz procedure [18]. In 
this procedure the evolution of the approximation error in m iteartions is represented by the multiplicative 
operator, P^u = I — Emu- Here, E^u = Km ' ' ' E^i^o ^^d Ki are coarse grid correction matrices. The 
classical domain decomposition approach does not work as a direct solver and therefore Pmu 7^ E Thus, 
the two-grid configuration in Fig. |7a] corresponds to a direct multiplicative Schwartz configuration. 

The following theorem establishes conditions on the inter-grid filters to obtain a direct solver. 

Theorem 4: A red-black harmonic two-grid configuration arranged as shown in Fig. |7al with non- 
singular A, A, Fji, A and Fj, works as a direct solver; i.e., v = u, if and only if 

njj,LALn7,L = nR^H^^Hfll,H ■ (43) 

Proof: From (l42l ). the direct solution is obtained if and only if 

KK=[0 0]. (44) 

Using Theorem [3l the conditions in the diagonal blocks, which guarantee a perfect recovery of the 
solution, give 

= nij,i^njj,,L(AA)-in/,Ln/,jyALAj^ and (45) 

= nR^L^R^H{^^)~^'^I,Ht^I,L^L^H ■ (46) 

The conditions in the off-diagonal blocks, which take care of aliasing cancellation, give 

(I - Ili^L^-^IlR^L^yL)t^I,L^~^t^R,H^H 

= {I- ni,H^-^'nR,H^H)fli,L^-^flR,H^H and (47) 

(/ - ni,H^~^IlR,HAH)ni,HA-^fiR,LAL 

= (I - ni,L^-^nR,L^L)fll,HA-'nR,LAL . (48) 

Here, all matrices are diagonal and therefore their products commute. None of the eigenvalues of inter-grid 
filters are zero because they are non-singular. Then, it is safe to multiply by any inverse of these matrices 
if necessary for simplifications. Since the coarse grid matrices are non-singular, each of the conditions 
above can be multiplied by AA. After this multiplication, using ( [36l ) and ( [37l ). the simplifications for the 
four equations independently give the same condition shown in ( |43l ). ■ 
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It is clear from (l43l) that inter-grid filters satisfying this condition should swap the low- and high- 
frequency eigenvalues of the system matrix. The mirror of the system matrix fits perfectly for this purpose. 
The following corollary gives a particular solution which tries to keep the algorithm as simple as possible. 

Corollary 2: A red-black harmonic two-grid configuration arranged as shown in Fig. |7al with non- 
singular A and such that det{AL + / 0, works as a direct solver with the following configuration: 

Fi = I, Fr = I, 

(49) 

Fi = A* and Fr = I . 

The coarse grid correction matrices for this configuration are 

K = {Al + Ah)-^ W 

. (51) 

Proof: Using the eigen-decomposition of filters and Theorem [2] gives 

fl/ = / , nR = I , 

(52) 

n/,L = Ah , Ui^H = Al and Ur = I . 

Using these eigenvalues in ( [36l ) and ( [371 ) gives A = + Ah and A = 2AlAh- Therefore, by Lemma 
[T] the coarse grid system matrices are invertible. Then, Theorem |4] can be applied and the eigenvalues of 
inter-grid filters above fulfill the condition ( |43] ). Finally, using ( [52l ) in Theorem [3] gives (|50l ) and dSTT ). ■ 

From (ISTI ) we see that the nested iteration step is not filtering nor reducing any frequency component 
of the error. It just equals the gain of all the effects. On the other hand, in (l50b we see that the correction 
scheme acts as a mirror filter by swapping the low- and high-frequency eigenvalues of the system matrix 
in the diagonal blocks. The aliasing effect in the off-diagonals is adjusted to cancel the symbols at the 
same row in K, so that KK = 0. 

This solution is particularly simple in the nested iteration step, since only down/up-sampling operations 
are used. The coarse grid matrix A = DAU has better sparseness than the system matrix A. The 
sparseness of A = DA*AU depends on the structure of down-sampling and non-zeros in A. 

When solving Poisson's equation, this solution is equivalent to the total reduction method by Schroder 
and Trottenberg [7], [8]. The derivation of total reduction methods was based on the structure of a 
stationary impulse response for LSI systems and thus imposes stronger assumptions on the system. 



Ah 
-Al 



-Ah 
Al 



and 



(50) 



K = -W 
2 



/ / 
/ / 
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(a) Multiplicative scheme 



(b) Additive scheme 



Fig. 7. Two-grid schemes approaches to solve Au — f exactly. The multiplicative scheme is different than the full two-grid 
scheme shown in Fig. [2a] because no smoothing iterations are used, and different coarse grids are used at each step. The additive 
scheme is the two-grid algorithm's version of the two-channel multi-rate system in Fig. [6] In both cases the two coarse grids 
form a partition of the fine grid and thus capture all the information from the fine grid. This fact allows these schemes to work 
as direct solvers; i.e., allows v = u. 



B. Additive Approach 

The second approach is the two-grid algorithm shown in Fig. |7b] which uses the structure of the two- 
channel multi-rate system from Fig. |6] This is a two-grid scheme with two nested iterations working 
in parallel at red and black coarse grids. The red and black coarse grids work as the two channels of a 
multi-rate system, but now, linear systems of equations are solved before the interpolation and addition 
of the two approximations. In a multi-rate system the idea is to reconstruct the input signal (in this case 
/) but now the idea is to transform this signal into the solution of the linear system. 

If this algorithm works as a direct solver then it means that the problem is being factorized into coarse 
grid sub-problems. The following proposition shows this in terms of a factorization of A~^. 

Proposition 4: The system shown in Fig. ITb] works as a direct solver; i.e., v = u, if and only if the 
following decomposition applies on the inverse of the system matrix: 

A-^ = IiA-^Ir + IiA-^Ir . (53) 

Proof: The approximation error in Fig. |7b]is given hy e = u — {v + v). Using ^ gives 

e= (^K + K - I^u (54) 

Then, an exact solution is obtained if and only if K + K = I. Using (l32l ) and (1331 ) gives (1531 ). ■ 
As opposed to (1411 ). here no cross-terms appear in the decomposition. This reflects the fact that nested 

iterations run independent of each other. This makes this scheme better suited for parallelization. 

This factorization is known in domain decompositions as the additive Schwartz procedure [18]. In 

this procedure the evolution of the approximation error in m iteartions is represented by the additive 
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operator, Pad = Km + ■ ■ ■ + Ki + Kq, where Ki are coarse grid correction matrices. The classical 
domain decomposition approach does not work as a direct solver and therefore Pad ^ I- Thus, the 
two-grid configuration in Fig. |7b] corresponds to a direct additive Schwartz configuration. 

The following theorem establishes conditions on the inter-grid filters to obtain a direct solver. 

Theorem 5: A red-black harmonic two-grid configuration arranged as shown in Fig. |7bJ with non- 
singular A, A, Fji, A and Fj, works as a direct solver; i.e., v = u, if and only if 

IlR,LflR,LAlUl,Lfll,L = UR,HffR,H^ltfl,Hfll,H , (55) 
UR,LffR,HAifi-I,Lfi-I,L + UR,HfiR,HAHfi-I,Hfil,L = 

IlR^HflR,LAlni,Lfll,L + nR^HflR,HAHtfl,Lfll,H and (56) 
nR^Lf^R,LAlt^I,Lf^I,H + nR^HflR,L^H^I,Hfil,H = 

tLR,LilR,LA.lni,Hlll,L + ILR,LllR,HAjiI[i,Hni,H ■ (57) 

Proof: From ( [54l ) the direct solution is obtained if and only if 

i? + K= [^0] . (58) 

Using Theorem [3l the conditions in the diagonal blocks, which guarantee a perfect recovery of the 
solution, give 

nj^L^-^'nR,L^^L + fli,L^'^flR,LAL = I and (59) 

nj^HA-^'nR,HAH + fll,H^~^flR,H^^H = I ■ (60) 

The conditions in the off-diagonal blocks, which take care of aliasing cancellation, give 

n/,LA-injj,H - n/.LA-^nH./f = and (6i) 

ni,HA-^tiR,L - Ui,HA-^flR,L = . (62) 

Here, all matrices are diagonal and therefore their products commute. None of the eigenvalues of inter- 
grid filters are zero because they are non-singular. Then, it is safe to multiply by any inverse of these 
matrices if necessary for simplifications. Since the coarse grid matrices are non-singular, each of the 
conditions above can be multiplied by AA. After this multiplication, using (l36l ) and (|37] ). the algebra on 
the conditions ( |59l ) and ( [60l ) simplifies to the same condition in dSSl ). and the algebra on ( [6l1 ) and ( [62l ) 
simplifies to ( [56l ) and dSTl) . respectively. ■ 



20 



IEEE TRANSACTIONS ON SIGNAL PROCESSING, PREPRINT, SEPTEMBER 2009 



These conditions are analogous to Vetterli's conditions (I23ti25l) . The analogy is more clear from ( [59l - 
[62l) where the only difference with (|23] - I25l ) is the existence of the matrices A, A, and Kh- This 
indicates the fact that a linear system of equations is being solved. 

Again, the mirror of the system matrix can be used to find a solution of these conditions. The following 
corollary gives a particular solution which tries to keep the algorithm as simple as possible. 

Corollary 3: A red-black harmonic two-grid configuration arranged as shown in Fig. |7bJ with non- 
singular A, works as a direct solver with the following configuration: 

Fi = A\ Fr = I, 

(63) 

Fi = A* and Fr = I . 
The coarse grid correction matrices for this configuration are 

and (64) 

y . (65) 

Proof: Using the eigen-decomposition of filters and Theorem |2] gives 
n/,L = Ah , fii,H = Al, flij = / , 

(66) 

n/,L = Ah , Ui^H = Al and Ur = I . 

Using these eigenvalues in (l36l ) and (l37l) gives A = A = 2AlAh- Therefore, by Lemma [U the coarse 
grid system matrices are invertible. Then, Theorem |5] can be applied and the eigenvalues of inter-grid 
filters fulfill conditions (15514571) . Using the eigenvalues from (l66l ) in Theorem [3] gives (l64l ) and (l65l ). ■ 

The coarse grid correction matrices in (l64l ) and (l65l ) equal the gain of filtering and aliasing effects. 
Interestingly, the symbols of K and K correspond to the black and red harmonic aliasing pattern in ( fT6l ). 
N and N, respectively. The aliasing effect in the off-diagonals have opposed signs between red and 
black coarse grids, so that K + K = I. 

Compared with the solution for the multiplicative approach, the additive approach involves more com- 
putations. This is because both red and black coarse grid matrices use a mirror filter. On the other hand, 
this approach is better suited for parallelization. Therefore, both the multiplicative an additive approaches 
become useful depending on the computational resources available. The multiplicative approach is more 
convenient in a single processor and the additive approach is more convenient in a multi-core architecture. 



K = -W 



I -I 
-I I 



K 



W 



I I 
I I 



NAVARRETE MICHELINI: DIRECT MULTI^GRID METHODS 



21 



VII. Direct Multi-grid Methods 

The convergence analysis of previous sections will be valid at each coarse level if harmonic aliasing 
patterns exist for each coarse system matrix. The following definition introduces the multi-scale property 
needed on the biorthogonal eigenvectors of the system. 

Definition 5 (Multi-grid Harmonic Basis): A multi-grid harmonic basis of level is any biorthogonal 
basis of C". A multi-grid harmonic basis of level / > is a biorthogonal basis of C", n = 2' no and 
no G N"*", with red-black harmonic aliasing patterns and such that the down-sampled low-frequency 
eigenvectors form a multi-grid harmonic basis of level I — 1. 

The two-grid methods in section |Vl] can be extended to multiple grids if the eigenvectors of a system 
matrix A G C"^", with n = 2'no and no = 0(1), form a multi-grid harmonic basis of level I. 

A. Multiplicative direct multi-grid algorithm 

The recursive implementation of the multiplicative approach shown in Fig. |7a] is explained in Table J] 
Here, the difference equation for the number of multiplications performed gives m(n) = 0{n\ogn). 

Starting from two coarse grids in the first level, the algorithm creates coarser grids which form a 
partition of the fine grid into more subsets of nodes. In Fig. [8a] the sequence of coarse grids visited by 
the algorithm is shown for / = 3. The structure is a W-cycle, well known in multi-grid methods [16]. At 
each level the partition of nodes is duplicated. For instance, in the coarsest level the grid partition gives 

where the subindexes from left (fine) to right (coarse) indicate if red or black grids were chosen. 

This algorithm has the nice property that the red coarse system matrix A = DAU reduces the sparseness 
of the system matrix A. In coarser levels the system matrix might soon become diagonal and the system 
is solved in linear time. Thus, there are good chances that the structure in Fig. |8a] changes from a W-cycle 
into a V-cycle [16], reducing the computational complexity from 0{nlogn) to 0{n). This is actually 
what happens when solving Poisson's equation, where this method is equivalent to total reduction methods 
[7], [8]. In general, this depends on the structure of the system and it will not be study here in depth. 
The example in section [Villi will show a case where the computational complexity is effectively reduced. 

B. Additive direct multi-grid algorithm 

The recursive implementation of the additive approach shown in Fig. |7b] is explained in Table Jl Here, 
the difference equation for the number of multiplications performed gives m{n) = 0{nlogn). 
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TABLE I 

Pseudocode for the direct multi-grid algorithm following multiplicative and additive approaches to 
SOLVE Au = f. The matrix A* represents the mirror of A according to definition[3J which only involves 
changing the sign of some of the entries in A. It follows from both cases that the number of 

multiplications IS m(n) = 0{n\ogn). 



Task 



Multiplications 



u — DMG_multiplicative(yl, /) m(n) 
if n > no , 710 — 0(1) 

f = Df 

A = DAU 

V = DMG_multiplicative(A, /) rn{^) 
vo — tJv 
r^f~Av 0{n) 
r = Dr 
Ii = A*U 
A = DAIi 0(n) 

= DMG_multiplicative(j4, r) '"(f) 

eo = Iiv C'('^) 

V = Vo + eo 
return(t;) 

else 

return(A-V) 0{1) 



Task 



Multiplications 



u = DMG_additive(A, /) 

if n > no , no = 0(1) 
f = Df 
h = A*U 
A = DAIi 



Ii = A*U 
A = DAIi 



U = IjV + IlV 

return(M) 
else 

return(A"\f) 



m(n) 




0(n) 



V = DMG_additive(A, /) m(f ) 





0(n) 



V = DMG_additive(y4, /) m(f ) 



0(n) 



o(i) 



Same as for the multiplicative approach, the algorithm creates coarser grids which form a partition of 
the fine grid into more subsets of nodes. In Fig. [8b] the sequence of coarse grids visited by the algorithm 
is shown for / = 3. The structure is a binary tree that duplicates the number of partitions at each level. 
The partition of the coarsest grid is the same as in (|67] ). Here, each of the coarsest grids is visited only 
once, which is possible because the solutions from different partitions do not depend on each other. 

This algorithm is more convenient for parallelization. In fact, from the finest to the coarsest grid in Fig. 
[8bl only down-sampling operations are performed. From the coarsest to the finest grid, only interpolations 
are performed. Therefore, these operations can be carried out in a single step and then the algorithm 
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(a) Multiplicative algorithm (b) Additive algorithm 



Fig. 8. Diagram of the sequence of computations in a direct multi-grid scheme following multiplicative and additive approaches. 
Moving from fine to coarse levels, each box adds an V or '6' at the right to indicate whether the red or black coarse grid 
is used, respectively. In the multiplicative approach, all the possible coarse grids are visited sequentially through a W-cycle in 
order to obtain the exact solution. In the additive approach, all the possible coarse grids are visited in parallel through a binary 
tree in order to obtain the exact solution. 




Fig. 9. Multi-grid scheme following a multi-channel additive approach to solve Au = /. Compared with the binary tree path 
shown in Fig. [Sb] here the source vector / is directly taken to the coarsest levels. The coarsest grids form a partition of the fine 
grid with the eight possible combinations of red (r) and black (b) coarsening. 

becomes analogous to a multi-rate system with multiple channels. The diagram of this algorithm is shown 
in Fig. |9] The computational complexity of a parallel implementation divides 0{nlogn) by the number 
of channels and reaches 0{logn) in the best scenario. 

VIII. Example 

Consider a two-dimensional system in which A is the finite-difference discretization of Helmholtz's 
equation — V^u — k'^u = on a square domain with periodic boundary conditions and unit step size. 
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(a) J7 = n,. + Qf, 



(c) ^ — ^frf ~\-^j-fb~\~^rbr ~\~^rbb~\~^brr ~\~^brb~\~^bbr ~\~^bbb 

Fig. 10. Red-black partition of nodes in a 2D square domain, Dark squares indicate the selected nodes. 
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(a) Impulse response in Q 
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(b) Impulse response in ilr- 



(c) Impulse response in fir 



Fig. 11. Impulse response coefficients from a finite-difference discretization of the Laplacian operator: — V^. The coarse grid 
system matrices obtained with mirror filters give impulse responses that grow in coarse grids. 



The size of the problem is set to n = 32 x 32 and the wavenumber is set to A; = ^. Thus, the stationary 

-1 

-1 4-fc2 _i (the underUne denotes the diagonal element). 



impulse response of A is 



-1 



The system is invariant under space shifts and therefore the eigenvectors are given by harmonic 



functions {W)ij = exp 



4 

rj, 



After proper reordering, the basis has harmonic aliasing patterns 



for the down-sampling shown in Fig. [TOl The chessboard down-sampling shown in Fig. llOal gives a 



mirror matrix A* with impulse response 



1 4-fc2 1 



An important problem arises naturally for systems in two or more dimensions. The down/up-sampling 
by a factor of 2 is not enough to reduce the impulse response of the product AA*. Then, the impulse 
response of system matrices grows larger and larger in coarse grids. An example for the impulse response 
of a Laplacian operator is shown in Fig. [TT] The same problem appears in other direct multi-grid solvers 
like total and partial (or cyclic) reduction methods [7]-[9]. 
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(a) Source vector: /. (b) Exact solution: u. (c) Source vector: /. (d) Exact solution: u. 



Fig. 12. Helmholtz's equation, ~V^u~ k^u — 0, with fc = ^, is solved in a square domain with periodic boundary conditions 
for two source vectors. In Fig. I12al the source vector is set to f{i,j) ~ sin(i7r/16) sin(j7r/16) + sin(i7r/2) siii(j7r/2), mixing 
low- and high-frequency sources. The exact solution is shown in Fig. I12bl In Fig. I12cl the source vector is equal to 1 in four 
neighboring nodes at the center of the figure and zero elsewhere. The exact solution is shown in ll2dl 




(a) vo for / in Fig. I12al (b) eo for / in Fig. I12al (c) vo for / in Fig. I12cl (d) eo for / in Fig. I12cl 

Fig. 13. Intermediate solutions of Helmholtz's equation using the direct two-grid algorithm with a multiplicative approach. 



In Fig. I12al and I12bl a solution of Helmlioltz's equation is shown when the source vector is the 
superposition of two frequencies, low and high. In Fig. I13al and Il3bl the intermediate solutions of a 
two-grid multiplicative approach are shown. In this case the red coarse system matrix is diagonal and the 
high-frequency in the forcing function is not seen at the red coarse grid. Therefore, the solution from the 
red coarse system is trivial and only considers the low-frequency component. The black coarse system 
adds a coiTcction with the high-frequency and part of the low-frequency component of the solution. In 
Fig. I14a[ I14cl and I14el the intermediate solutions of an additive multi-grid approach are shown for three 
coarse levels. The high-frequency component is not seen in some of the coarse grids. 

In Fig. I12cl and I12dl a solution of Helmholtz's equation is shown for a sparse source vector. In Fig. 
|13c| and |13d| the intermediate solutions of a two-grid multiplicative approach are shown. Since the red 
coarse system matrix is diagonal, the solution of nested iteration is trivial, but not very significant as the 
correction scheme gives most of the solution. In Fig. I14b[ I14dl and I14fl the solutions at coarse grids are 
shown by using an additive approach. One of the non-zero values of the source vector appears in each 
of the coarse grids ^rr, ^rb^ and 0;,;,. At the third coarse level the source vector is zero in four of 
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(e) Wrrr, Urrft, Urbr, Urbb, Ubrr, Ubrb, Ubbr and Ubbb 



(f) Urrr-, f^rrb^ '^rbrt '^rbh-! T-^brrt '^brhi Ufjhr 

and Ubbb 



Fig. 14. Intermediate solutions of Helmholtz's equation using the direct multi-grid algorithm with an additive approach. 
Considering the two sources in Fig. |12a| and |12c| the exact solution is given by: u — Ur + Ub al the first coarse level in 
Fig. |14a| and |14b| respectively; u — Urr + Urb + Ubr + Ubb at the second coarse level in Fig. |14c| and |14d| respectively; and 
u = Urrr + Wrrb + Urbr + Urbb + Ubrr + Ubrb + Ubbr + Ubbb at the third coai'se Icvel in Fig. 1 14el and 1 14f I respectively. 



the eight grids. This is a consequence of the down-sampling interpolation and shows the advantage of 
this approach when the source vector is sparse. 

IX. Conclusions 

Numerical methods to solve linear systems of equations were obtained based on the similarities of 
the full two-grid algorithm and perfect reconstruction filter banks. The two alternatives, multiplicative 
and additive, correspond to direct Schwartz domain decomposition methods based on a partition of the 
original domain. The additive approach can be used to parallelize the problem among all the available 
processors, whereas the multiplicative approach is more efficient in a single processor. 

Future research will focus on the application of these algorithms in systems that are not LSI. On one 
hand, the algorithms are ready to work on systems that are known to have harmonic aliasing patterns 
but more numerical studies are necessary. And, on the other hand, the most challenging problem is to 
understand the physical and geometrical implications of harmonic aliasing patterns. This is essential 
to construct practical methods to check aliasing patterns and be able to use these methods in more 
challenging problems. 
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Appendix A 
Proof of Theorem [H 

Proof: The proof that the red harmonic aUasing pattern is equivalent to ( fTSl ) is due to Navarrete and 
Coyle [11]. The proof that the black harmonic aliasing pattern is equivalent to ( fT9l ) follows the same 
reasoning. Given the partition of eigenvectors W = \WlWh\ and V = [VlV^], the black harmonic 
aliasing pattern can be written as the following set of biorthogonal relationships 

{DVLfiDWL) = i / , (68) 

{DVl)''{DWh) = -\I , (69) 

{bVH)"{DWL) = -i / and (70) 

{DVh)''{DWh) = \I . (71) 

Since W and V form a biorthogonal basis, WV^ = WlV^ + WhVh = ^- Pre-multiplying by D and 
post-multiplying by U gives 

{bw){DV)^ = {bwL){bvL)" + {bwH){bvH)^ = i . (72) 

First, ([T9l ) is assumed. Then, equation ( 1721 ) immediately implies the set of biorthogonal relationships 
above, and the black harmonic aliasing pattern is fulfilled. Second, the black harmonic aliasing pattern 
is assumed. Pre-multiplying ( 1721 ) by {DVl)^ and using equations ( |68l ) and (|69l ) gives DVl = —DVh- 
Similarly, post-multiplying ( 1721 ) by DWh and using equations ( |69l ) and ( TtTI ) gives = —DWh- 

Therefore, the black harmonic aliasing pattern implies (|T9l ). ■ 

Appendix B 
Proof of Lemma [H 

Proof: The proof of (l34b is due to Navarrete and Coyle [11]. The proof of (l35l) follows from 

I-i = |i3FR^F7f7}~^ (73) 

= {{DW)URAUj{DVfy^ (74) 

= {(Wl)A(5Fl)^}"' (75) 

= 4. {bWL)A-\DVL)" , (76) 

where ( |30l ) is used in (1731 ). the eigen-decompositions of filters is used in ( 1741 ). Theorem [T] is used in (1751 ) 
and the biorthogonal relationships ( |68l ) to (TtTI ) are used in (1761 ). ■ 
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Appendix C 



Proof of Theorem [3] 



Proof: The proof of (l38l) is due to Navarrete and Coyle [11]. The proof of ( [39l ) follows from 



K = I- FiUA-^DFrA 



(77) 



= / - 4 WUi{V"UDWl)A-^{V^UDW)UrAV" 
= WV^ - 4 T^fl/ (i [ ] ) (i [ / -/ ]) UrAV 



H 



(78) 



(79) 



1^ ■f-n/.iA-inH^i.Ai n7,tA-inj^,HAjj 



(80) 



where (l30l ) is used in (1771) . the eigen-decomposition of filters and Lemma [T] are used in (1781 ). definition 
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